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' ~~* I Abstract. We examine the bar instability in models with an exponential disk and a cuspy 

I ""^ . NFW-like dark matter (DM) halo inspired by cosmological simulations. Bar evolution is studied 

' as a function of numerical resolution in a sequence of models spanning 10*~^ DM particles - 

^•J [ including a multi-mass model with an effective resolution of 10^" . The goal is to find convergence 

^^1 in dynamical behaviour. We characterize the bar growth, the buckling instability, pattern speed 

1^^ ' decay through resonant transfer of angular momentum, and possible destruction of the DM halo 

r^ , cusp. Overall, most characteristics converge in behaviour for halos containing more than 10^ 

CIh' particles in detail. Notably, the formation of the bar does not destroy the density cusp in this 

i ' case. These higher resolution simulations clearly illustrate the importance of discrete resonances 

t, , in transporting angular luomentum from the bar to the halo. 

t/j . Keywords, galaxies: spiral, structure, evolution; dark matter, methods: n-body simulations. 
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1. Introduction — The bar instability in a cold gravitating disk plays a major role in 
a spiral galaxy's dynamical evolution. At least 2/3 of spiral galaxies host bars (Knappen 
et al. 2000) and the fraction has not evolved significantly since z ~ 1 (Jogee et al. 2004; 
Sheth et al. 2008) . As models of galaxy formation become more sophisticated and reveal 
complex dynamical behaviour, it is important to understand the details of different phys- 
f^ ical processes that shape their morphology as well as to verify that numerical resolution 

OO is in fact adequate to follow their evolution. The bar-halo interaction is the driving mech- 

^^ • anism in disk galaxy evolution. As a bar churns through the DM halo with a pattern 
^ , speed 51b resonant interactions with halo orbits - a form of dynamical friction - transfer 
angular momentum from the bar to the halo and cause it to spin down (Tremaine & 
Weinberg 1984). This process was first pointed out by Lynden-Bell & Kalnajs (1972) 
C^ ' and has been studied in models with idealized rigid bars (Weinberg 1985; Hernquist & 
Weinberg 1992; Weinberg & Katz 2002; Weinberg & Katz 2007) as well as in models 
in which a stellar bar forms self-consistently in an unstable disk (e.g., Sellwood 1980; 
Debattista & Sellwood 1998; O'Neill & Dubinski 2003; HoUey-Bockelmann et al. 2005; 
Martinez- Valpuesta et al. 2006). There has been some concern that the process is too 
efficient, leading to bars that are much smaller than their corotation radii and so dis- 
crepant with observed bar galaxies (Debattista & Sellwood 2000). More recent studies 
with greater resolution suggest that the bars tend to lengthen moving out to their co- 
rotation radii as they slow down and so perhaps they are not inconsistent with reality 
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(O'Neill & Dubinski 2003; Martinez- Valpuesta et al. 2006). Weinberg & Katz (2002) hold 
a cautious view that lower resolution simulations can lead to spurious results because of 
the diffusive nature of noise that may move orbits into and out of resonances artificially 
while insufficient particle numbers may also underpopulate the resonant regions of phase 
space. While most current work to date has used ~ 10^ particles, they claim that as 
many as 10^ DM particles (or more) may be needed to simulate the resonant process 
with iV-body methods. 

In this study, we attempt to clear up the inconsistencies of current work and address 
the problem of numerical resolution hoping to converge to the correct physical behaviour. 
We present a series of bar-unstable disk-|-halo A^-body models with increasing resolution 
spanning a range Nh = lO** — 10* DM particles with Nd = 1.8 x 10^~^ disk particles. 
One further simulation uses a multi-mass method that increases the halo particle num- 
ber density by 200x in the halo centre so giving an effective Nh ^ 10^^°. The mass 
model is constructed from a 3-integral distribution function (Widrow & Dubinski 2005) 
describing an exponential disk embedded within a tidally truncated NFW halo and is 
similar to the model studied by Martinez- Valpuesta et al. (2006). Natural units for the 
simulations are Z? = 10 kpc, M = 10" Mq, F = 208 km/s, and T = 47.2 Myr. We 
discuss results both in simulation and physical units throughout and clarify when nec- 
essary. (For further details on the models and simulations see Dubinski, Berentzen & 
Shlosman 2008). Animations of the simulations are available for viewing at the URL: 
www.cita.utoronto.ca/~dubinski/IAU254 along with higher-resolution figures. 

Animation 1 shows the evolution of the disks in six models with increasing resolution in 
face-on and edge-on views. The lowest resolution simulations with Nh ^ 10^ are clearly 
deficient and either lose the bar or suffer from heating affects. At higher resolution, 
the behaviour is similar exhibiting the buckling instability and relaxation to a bar in 
quasi-equilibrium that gradually lengthens and slows down. Animations 2 & 3 show the 
multi-mass model with Nh = 10* in an inertial and corotating frame and illustrate how 
the bar grows from noise from the inside out saturating as a thin bar on reaching the 
corotation radius and then evolving into a fatter bar with a peanut-shaped bulge after 
the buckling instability. 



2. Bar Grovifth and Pattern Speed Evolution We study bar growth using the 
normalized Fourier amplitude \A2\ of the m — 2 disturbance in the plane of the disk 
within a fixed radius R < 0.5 (5 kpc). Figure [T] shown the exponential growth of 1^2 
before saturation. Higher resolution simulations reach saturation at later times. This time 
delay is the result of the lower amplitude of the Poisson fluctuations that seed the bar. 
Since the instability grows from these fluctuations if takes longer for them to saturate 
if the initial amplitude is lower. We estimate the time delays from the peak in \A2\ and 
synchronize the simulations for comparison of various evolutionary characteristics. 

We can also use A2 to measure the phase angle of the m = 2 mode and so estimate the 
pattern speed by taking the difference between subsequent snapshots. Figure [2] shows the 
pattern speed evolution as a function of numerical resolution. There is some scatter in 
behaviour that can be accounted for from the expected variance introduced by different 
initial conditions but overall the agreement is good. The highest resolution simulations 
show a modulation of the pattern speed at a frequency close to fit itself. This probably 
indicates an interaction between the bar and the gravitational wake in the halo that only 
shows up with sufficient numerical resolution. 
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Figure 1. Initial growth of the m — 2 Fourier component \A2\ for stars with R < 0.5 for 
two model sequences using Nd = ISA", 180K, 1.8M, 18M with Nh = WOK, IM, lOM, lOOM re- 
spectively. The In jyl2| grows approximately linearly with time independent of the choice of Nd 
and Nh showing the exponential growth of the bar mode. The dashed line shows an exponen- 
tial timescale that is approximately r = 8 (370 Myr). Since the bar grows from the Poisson 
noise within the disk then we expect the noise amplitude to be proportional to N^'^ so larger 
simulations will saturate at later times. 

Figure 2. Evolution of the pattern speed Q,b for two model sequences at different resolution. 
The curves have been shifted in time so that the bar growth evolution is coincident with the 
IM particle model. 



3. Halo Density Profile — We also measured the evolution of the halo density profile 
over the course of the simulation. Previous work has shown both preservation and destruc- 
tion of the density cusp so we focus on the processes in the central regions. Only studies 
using the self-consistent field A^-body method (SCF) lead to a core (HoUey-Bockelmann 
et al. 2005; Weinberg & Katz 2007b) and there has been some concern about numeri- 
cal instabilities that arise in those methods (Selwood 2003). Figure [3] exhibits the final 
density profile along with the difference from the initial profile as a function of halo par- 
ticle number. The profiles agree well to within the limit of their gravitational softening 
radius and show convergent behaviour. The cusp is not destroyed in this case but rather 
the central density increases modestly, by a factor of 2, as a result of the bar evolution, 
following its buckling, which leads to an increase in the central stellar density (Sellwood 
2003). 

The bar that forms in this simulation has a similar mass ratio Mt/Mhaio ~ 0.6 but 
is much thicker than the fiducial rigid bar simulation in Weinberg et al. (2007b) that 
destroys the cusp. Their study also included thick bars which did not destroy the cusps 
and the bar that forms here overlaps with those in their study. We conclude that there is 
no direct contradiction with the most recent results of Weinberg et al. (2007b) but would 
argue that the thicker bar models are probably more relevant to real galaxies. Thin bars 
do not persist for long before responding to the buckling instability and so the rigid bar 
approximation is not applicable over a Hubble time and probably not relevant to most 
galaxies (Martinez- Valpuesta & Shlosman 2004). 



4. Orbital Resonances — The bar slowdown is the result of dynamical friction that 
leads to angular momentum transport to the DM halo. The process is due to resonant 
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Figure 3. A comparison of density profiles at f = 7.08 Gyr for different halo particle numbers 
Nh- We also show the initial density profile (dashed line) and the best fit NFW model curve 
(dotted line) to the initial profile over the range < r < 100 kpc. The NFW parameters for 
the fit are rs = 4.3 kpc, Vmax ~ 160 km s~^, where Vmax is the maximum circular velocity 
at r = 2.16rs = 9.3 kpc. Note that this halo is more concentrated than the typical galactic 
dark matter halos in cosmological simulations to model the contraction expected during dissi- 
pative galaxy formation. The dotted vertical lines show the softening length e used at different 
resolutions. 



interactions between the rotating bar's pattern speed and the halo particles' azimuthal 
and radial orbital frequencies. When hflr+h^rp = mflb then orbital resonances occur and 
halo particles torque or are torqued by the bar leading to a change in angular momentum. 
We can estimate the orbital frequencies at a specific time by freezing the potential and 
integrating the orbits of particles in this potential rotating with the pattern speed at that 
time (e.g., Athanassoula 2002; Martinez- Valpuesta et al. 2006). Spectral analysis can then 
be applied to the orbital time series to determine the fundamental orbital frequencies 
r^r and Q,ij, (Binney & Spergel 1982). (Note we label these frequencies with the usual 
epicyclic variables k = fJ^ and fJ = U,^ in our figures below.) The dimensionless frequency 
rj = {Q.~Vti,)/ K, is a useful way to characterize resonances since the values 77 — 1/2, 0, —1/2 
correspond to the inner Lindblad, corotation, and outer Lindblad resonances for m = 2. 
Further negative half-integer values correspond to higher order resonances that can also 
absorb the angular momentum. 

Figure[3]shows the change in the z component of the angular momentum J^ for particles 
binned as a function 77 between t — 100 and i = 150 at different numerical resolution. 
The spikes at the half-integer values reveal the resonances. Most angular momentum is 
transferred through the corotation resonance though it is also appears to be transferred at 
other frequencies, but more randomly — this difference is probably the result of particles 
in resonance before or after t = 100. The detailed behaviour of the distributions seem to 
converge for N^ ^ 10^. At lower resolutions, the resonant spikes have a smaller amplitude 
and higher order resonant interactions are missing. 
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Figure 4. The distribution in tiie ciiange in the z angular momentum AJz between t — 100 (4.7 
Gyr) and t — 150 (7 Gyr) plotted as a function of the dimensionless frequency 77 = (n — Q.b)/K, 
measured at f = 100. The spikes at r; = 0.5, 0.0, —0.5 correspond to the ILR, COR, and OLR 
respectively while other spikes refer to higher order resonances. This plot shows how halo par- 
ticles in resonant orbits are the main sink of angular momentum. The detailed distributions are 
converging for Nu > 10^ particles while lower resolution simulation miss some of the higher 
order resonances. 



Finally, we examine the change in halo phase space density by computing the particle 
number density in {E, L^) space and computing the difference between i = and t = 150 
in a similar way to HoUey-Bockelmann et al. (2005). In this way, we clearly see the 
resonant regions visible as discrete islands of particle overdensity in (S, L^) space (Fig. [5]). 
We can also overplot the values of {E, L^) for the particles found in the resonant spikes 
in the analysis to see where they lie in phase space. The right panel of Figure [S] clearly 
shows that these islands are directly related to the discrete resonances extracted in our 
spectral analysis. Animation 4 describes the time evolution of the differential number 
density in phase space and reveals how the resonances move through a large fraction of 
the halo mass. By counting particles in resonant peaks at different times we estimate 
that roughly 30% of the halo particles are affected. 

We conclude that the resonances are broader than thought and so simulations with 
more than IM halo particles do a reasonable job of tracking bar evolution. However, a look 
at the distribution of orbital frequencies reveals that higher order resonances are missed 
at lower resolution with less than lOM particles. This effect could account for the different 
rate of angular momentum loss at higher resolution. Future studies should examine the 
bar instability self-consistently using the same initial conditions with different iV-body 
methods to resolve current inconsistent results on the cusp/core evolution of DM halos. 
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Figure 5. Change in particle number density in {E,Lz) space between t = and t = 150 (7.0 
Gyr) for the A^^ = 10* single mass model. The resonant regions show up clearly as peaks (red 
regions) in phase space in the left panel. In the right panel, we overplot the {E, Lz) coordinates 
of a random subset of particles located at discrete resonances at i = 150 within Srj = ±0.05 
(black-ILR-7) = 0.5, red-COR-77 = 0.0, green-0LR-?7 = —0.5, blue-77 — —1.0, magenta-77 — —1.5, 
and cyan-r; = — 2.0 
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